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Abstract 

The principle of stationary phase (PSP) is re-examined in the context of linear time-frequency (TF) decomposition 
using Gaussian, gammatone and gammachirp filters at uniform, logarithmic and cochlear spacings in frequency. This 
necessitates consideration of the use the PSP on non-asymptotic integrals and leads to the introduction of a test for 
phase rate dominance. Regions of the TF plane that pass the test and don't contain stationary phase points contribute 
little or nothing to the final output. Analysis values that lie in these regions can thus be set to zero, i.e. sparsity. 
In regions of the TF plane that fail the test or are in the vicinity of stationary phase points, synthesis is performed 
in the usual way. A new interpretation of the location parameters associated with the synthesis filters leads to: (i) 
a new method for locating stationary phase points in the TF plane; (ii) a test for phase rate dominance in that 
plane. Together this is a TF stationary phase approximation (TFSFA) for both analysis and synthesis. The stationary 
phase regions of several elementary signals are identified theoretically and examples of reconstruction given. An 
analysis of the TF phase rate characteristics for the case of two simultaneous tones predicts and quantifies a form 
of simultaneous masking similar to that which characterizes the auditory system. 

Index Terms 

method of re-assignment; cochlear filters; gammatone; gammachirp; simultaneous masking; 



I. Introduction 

The principle (or method) of stationary phase (PSP) [ 1 ] is a result from asymptotics that can provide closed- form 
approximations, in the limit as A — > oo, to often intractable oscillatory integrals of the form 

/oo 
a(t)e> Xb ®dt (1) 
-oo 

where a, b,t, X £ IBL There are two steps involved: (i) recognition that in the limit the integral will be almost 
zero everywhere in the interval {t : — oo < t < oo} except near values of t where the derivative b(t) is zero, the 
stationary phase points; (ii) the integral in the vicinity of these stationary phase points can be expressed in terms of 
the second derivative of the phase i.e. b(t). Perhaps the most successful application of the PSP in signal processing 
has been in the context of synthetic aperture radar (SAR), where it is the starting point in the development of many 
of the Fourier-based imaging algorithms, c.f. O. Application of the PSP in not without its pitfalls. It is tempting 
to use the PSP in the non-asymptotic cases where A = 1 to find closed form approximations to integrals such 
as Fourier transforms. The argument for this requires that the phase b(t) is changing much more rapidly than the 
amplitude a(t). However, as pointed out in 0, a degree of care must be exercised, particularly with step (ii). 

The primary interest here is its application to linear time-frequency (TF) decomposition [4J. The motivation is the 
recent resurgence of interest in analogue filter banks both as part of a synthetic cochlea and as a means to provide 
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power efficient implementations of analysis filter banks |5J. The desire with both is to extract salient features from 
the TF decomposition using the limited functionality associated with analogue circuitry. This does not deny the 
considerable work that have been done on computational modelling of the auditory system, typified by papers such 
as (6j and the references therein. However the main purpose of such work is to model and predict the response of 
the auditory system to stimulus rather than to expose the signal processing principles that might be at work. 

The PSP is a natural place to start because of the prevalence of oscillatory terms in TF decompositions and the 
hope that the stationary phase points may provide a means for identifying salient features as well as a focus for 
sparse decompositions without the need for the usual iterative re-synthesize steps [7]. The PSP has been applied 
to linear TF decomposition for both analysis, (H & O, and synthesis ifTUll . the latter leading to the method of 
reassignment. The approach here is to revisit ifTOIl and to fundamentally re-interpret it to provide a PSP-based 
approximation to the TF synthesis integral. There is no attempt to either reassign iflOl or relocate |[TT1l components 
in the TF plane because of the limited functionality mentioned above. Another alternative would be to follow an 
amplitude-based approach such as ridgelets [12]. However this leads to algorithms that are far from feasible with 
analogue circuitry and further, ridges in the TF plane may not be appropriate when dealing with auditory filters 
such as the gammatone [13] and gammachirp lTT4l which have asymmetrical impulse responses and, in the case of 
the latter, an asymmetrical frequency response. 

Concentration on the synthesis rather than the analysis integral is advocated because: (i) most methods for sparse 
atomic compensation, c.f. 0, are based on re-synthesis of the original waveform; (ii) one of the main functions 
of the auditory system is to code the incident waveforms and coding requires at least some consideration of the 
potential for reconstruction (even when reconstruction is not a requirement); (iii) smooth variations of the magnitude 
and phase of the integrand are more readily satisfied for the synthesis integral than the analysis integral (because 
the signal has already been filtered in the analysis process). 

In Section [n] a preliminary description of the filter banks to be considered is given. A framework is presented 
that accommodates Gaussian, gammatone and gammachirp filters at uniform, logarithmic and cochlear spacings in 
frequency. More details of the relevant properties are provided in Appendix [A] Section III is an examination of 



the application of the PSP to non-asymptotic integrals and proposes the use of step (i), above, without step (ii). 
The concept of phase -rate dominance is introduced to partition the interval into sets where the PSP is and is not 



applied. The synthesis double-integral is addressed in section IV and the PSP is applied to it in a new way to 
facilitate the detection of stationary phase points in the TF plane. Extending the single-integral concepts to double 
integrals (Appendix IB]) leads to a test for phase-rate dominance in the TF plane. The stationary phase regions of 



several elementary signals are identified in section [V] and examples of reconstruction provided. Section VI examines 
the phase rate characteristics in the TF plane for the case of two simultaneous tones . It is shown that if these 
characteristics are used to test for the presence of a tone at a particular frequency, a form of simultaneous masking, 
similar to that which characterizes the auditory system [15], is observed. All of the above is addressed from the 
perspective of deterministic signals. Consideration of random processes is left for future work. The effects of noise 
on the related method of reassignment are covered in [16] and] 17]. Finally, in section VII, conclusions are drawn. 

II. Preliminaries 

Consider a time-frequency analysis X^it) of a signal of interest x(t) of the form: 

X u (t) = x(t) * hu(t), {UJ : Umin < U < Umax} (2) 

where * denotes convolution and the impulse response of a single filter in the analysis filter bank is given by: 

Kit) = ph(fit)&* (3) 

Each filter is formed using a prototype filter h(t) and a nominal bandwidth f3. The frequency response of the 
analysis filter is related to the frequency response of the prototype in a straight forward way, i.e. H u (fl) = 

F[hu(t)] = H (^r)- where H(Q) = F[h(t)\ = /f^ ^(*)e~ jm dt. Prototype filters that will be considered here 
are a gammatone pulse, a gammachirp pulse and a Gaussian pulse. While the Gaussian pulse is a common 
and analytically convenient choice for time-frequency analysis, the gammatone lfT3l and gammachirp pulses 11141 
more closely model the cochlea in the ear. Details of these prototype filters and their properties are summarized 
in Appendix [A] The prototype filters are normalised such that H (0) = 1. This property and multiplication by the 
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nominal bandwidth (3 in ^ normalizes the maximum gain of each filter to unity at u rad/s. A re-synthesis, x(t), 
of the signal of interest is performed using filters matched to h u (t). Thus for real signals of interest 

•i 



x(t) 



1 



Z(t, fj) drd/i 



(4) 



C wo 

with integrand Z(r,[i) = X^(r)/i* (r — t) and where C is a constant and and denote 'real part of 
and 'imaginary part of respectively. As a convenience, in order to deal with a number of possible filter bank 
spacings, define a filter bank variable fx that lies in the range [0, 1], A value \i = indicates the lower edge of 
the filter bank and fi = 1 indicates the upper edge. The frequency to at which the filters gain is a maximum is a 
function of fi and the bandwidth of the filter is proportional to the derivative i.e. f3 oc Hence the total 
number of filters required to just cover the band of interest is given approximately by In the following u m in 
is nominally the lowest frequency covered by the filter bank and io rnax is the maximum. For a uniformly spaced 
filter bank u = {u) max — Lu m i n }fj, + U) m i n and hence u = {u> — u> m i n }/ {^max ~ w m j n }, ^ = ui max — u) m i n , the 
nominal bandwidth (3 is a constant and 4^ = 0. For logarithmically spaced filter banks, similar to wavelets HI . 
^min^ b ^ where b = ln( ^ max ) and hence fi = ^ln ^, = few, the nominal bandwidth /3 is proportional 
is a constant. For a cochlear spaced filter banks, based on the approximation of JT8I for low sound 



to 

to to and 



dui 



pressure levels, u 

duj 



d/3 



+ a}e b ^ — a, where a 



2ttx10 3 



In 



and hence fi 



|1» 



qj+q 



4.37 ' u *~ \ ui min +a / 

j— — 6{w + a} and ^ is a constant. For all the above filter banks, the equivalent rectangular bandwidth (ERB) 
lfl31 and the 3 dB bandwidth are related to the nominal bandwidth (3 in a straightforward way. 



III. The PSP and non- asymptotic integrals 

Consider the integral |T]) evaluated over an interval [t— A/2, t+A/2]. Assuming that the function f(t) 
is well approximated by its Taylor series over this interval, gives 

/•t+A/2 

lA(t) ± / f(t')dt' 



a(t)e 



Jt-A/2 
fit) 



t+A/2 
t-A/2 



l i + /(*) 



l i + /(t) 2 



}di' 



The derivatives can expressed as 



and: 



/(*) 
/(*) 



/(*) 



a(t) ( a(t) 



d(t) 
a(t) 



+ j\b(t) 



a(t) \a(t) 



+ j\b(t) + 



d(i) 
o(t) 



+ jA6(i) 



In the asymptotic case, at a suitably large value of A, 



,/(*) 



j\b(t) and |g 



X 2 b 2 (t) + j\b(t) and hence the 



integral is dominated by the phase derivatives b(t) and b(t). Thus, at a stationary point, where b(t) = 0, the integral 
can be evaluated in terms of b(t) without reference to derivatives of a(t). For the non-asymptotic case where A = 1, 
these the approximations are also dependent on the relationships between the derivatives of the magnitude and phase 



of the integrand. Thus 



m 
tit) 



jb(t) provided, as in (8), that 

\b(t) 



d(t) 



o(t) 



(5) 



and f® 
and m 



-b 2 (t) + jb(t) provided also that 



6(t) 



> 



ii(t) 








fa(t)\ 


a{t) 


I a(t) J 




- 
dt 





(6) 
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If the amplitude and phase derivatives are available then at every value of t it is possible to test for what might 
be called first order or second order phase-rate dominance using ([5]) or, Q and ([6]), respectively. The approach 
adopted here, where a closed form approximation to the integral is not required, is to circumvent the difficulties 
associated with step (ii) of the PSP by simply not using that step. Rather Q is used to identify a set So that is the 
union of intervals of t where step (i) is valid and the compliment to that set So where it is not. If the stationary 
phase points {U}i where b(ti) = and intervals such as Si G {t : — 5\ < t < ti + 62,61 > 0,62 > 0} that 
contain them can be identified, the integral over the whole real line ([T} can be replaced by a similar integral over 
the union of intervals S G {{J { Si} U So- It is also worth noting that for ^ to be satisfied at or near a stationary 
phase point, would also require that |^y| — > as t — > to- Thus the normalized amplitude rate must go to zero 
more rapidly than the phase rate, i.e. ^ must apply. Thus there are liable to be be intervals where Si G So and 
for these intervals there is no need to identify 61 and 62 . 



IV. A Time-Frequency Stationary Phase Approximation 

The PSP can be extended to double integrals such as Q with the result defined in terms of the gradient and 
Hessian of the phase of the integrand, c.f. CQ, Pg- 478. However, as in Section [III] the full form the PSP is not 



used here. Specifically, the gradient is used in step (i) to identify the stationary phase points but the Hessian of 
step (ii) is not used to form an approximation to the integral in the vicinity of these points. A test similar to ([5]) is 
developed for double integrals in Appendix [B] to identify regions of phase-rate dominance in the TF plane where 
the PSP is applied by numerical integration in the vicinity of the stationary phase points. As in Section III it is 

convenient to define normalized time- and frequency- derivatives of the integrand Z(r,fi), i.e. Z T = s |^~y an d 



Z ^ = a "^( T '^ dM respectively, where 



z T = °\ r ;y + 



x u (t) K(T-t) 
&X u (t) h*(P{r-t}) 



+ /3ttt^ ttt (7) 



and 



where 



X u (t) ^h*((3{r-t}) 





t}\-j{r-t} 

(9) 

Both Z T and Z^ are additions of a signal dependent term, e.g. JiJ ^jyy-, and a filter dependent term, e.g. Jt ^rj^—j^~- 
The former is a function of the pair (ui, r) whereas the latter is a function of the pair (u, r — t). Note that (qj, r — t) 
are themselves parameters of the filter, specifically, the frequency where the filter has maximum gain u and the 
delay r — t between the input and output of the filter. In contrast to the method of re-assignment derived in iflOll . 
the delay term {r — t} is interpreted as the group delay 

{r-t}=g(u) 4 __d_Z^(0)| n=w (10) 

of the filter (|3]) at uj, where the notation ZH indicates the argument of the complex variable H. Simple expressions 
for the group delay are given in the Appendix [A] in terms of the group delays of the prototype filters at zero 
frequency. The justification for the use of ( 10 1 in Q, ^ and ^ proceeds as follows: {r — t} is the delay between 



the input signal x(t) and the output of the analysis filter X u (t); since each filter is tuned to have a maximum gain 
a particular frequency u, the delay through the filter is also tuned to the rate of change of the phase response at 
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the frequency where the gain is maximum. A particular filter is thus jointly labeled with both its frequency oj and 
the group delay at that frequency g{oj). Thus (jT]) and ^ can be written as 



few 

XJt) 



+ f3ri-juj (11) 



and 



z. 



&X u (t) , d/3 [1 , . , J __,. A l da; 



I X w (r) da; [ /3 J J d/i 

(12) 

respectively, where 77(0;) = |4 frffe]T T- Then because ^{^(a;)} = for all three filter types including the complex 
gammachirp (c.f. Appendix Hj, the time derivative of the phase of the integrand is 

^ - <13) 

and the frequency derivative is 




Time and frequency derivatives of the analysis integral (|2|) are constructed using the derivative filters -^.h^{r) and 
^K{t) respectively, c.f. HI, API : 

J^(r) = x(t) * ^K{t) (15) 

— X u (t) = x(r)*— ^(r) (16) 
Alternatively, as suggested in ifTOl . they can be derived from the output of ([2]) by direct differentiation of the signal 



X u (r) to obtain ( [15] ) and by using neighbouring analysis filters to obtain an approximation to (16 1. Stationary 
phase points {(ct>i,Tj)}j are solutions to: 

5{Z T } = 9{Z M } = (17) 

The derivative filters can be expressed in terms of the prototype filter h(t) and its derivative h(t ). E xpressions for 
the derivative h(t) of the Gaussian and gammachirp filters can be found in Appendix |a| From (17) there are two 



conditions that must be satisfied simultaneously for a stationary phase point to occur at (wj,Tj), specifically: 

1) the frequency, 8Z ^^ = 9 j ^-j^ j , observed at the output of filter at time r, is equal to the centre 
frequency oj of the filter; 

2) the delay, 9Z ^^ = 9 { a j^^ 1» observed at the output of the filter at frequency oj, is equal to the group 
delay g{oj) of the filter at that frequency. 

Together these define a signal matching condition: at (oj, r) the frequency and delay observed at the output of the 
filter must match the designed centre frequency and group delay of that filter. 

Locating stationary phase point requires a grid search over oj for a bank of analogue filters or over both oj and 
r, for a discrete-time filter bank. Such a grid search is not onerous since it is implicit in the implementation of the 



analysis integral. With a grid search there is always the risk of missing the pair (o^,^) that satisfy (17 1, This risk 
can be reduced by: (i) defining a phase gradient vector 

0(t,m) = MZ T } 3f{^}] T (18) 
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where the superscript T indicates matrix transpose; (ii) using the Euclidean norm of this vector to construct a test 
for stationary phase points, i.e. 

H(r,fx)\\ < Ci (19) 

where the threshold C\ is a small positive real constant. The Euclidean norm is used here for analytic convenience. 
Other vector norms may be appropriate and may have desirable properties with respect to ease of implementation. 

In addition to finding stationary phase points, equations ( [TTj ) and (12) can also be used to test for phase-rate 
dominance in the TF plane. For phase-rate dominance the inequality 

/Hi > ^ ™ 

N T > m)II 

must be satisfied, where the amplitude gradient vector is 

a(r,A*) = mZ T } nZ»}] T (21) 
and the threshold C2 is a positive real constant greater than or equal to one. The projection term 

ll a (^)ll 

is formed from the sum of the projection of the phase rate vector in the direction of the normalised amplitude rate, 
i.e. T a/||a||, plus the projection in the orthogonal direction, i.e. cf) T J] a/||a|| and hence W = [} l }]. The 
test is derived in Appendix |B| Thus the PSP divides the time frequency plane into two regions: a region So where 



( 20 1 is satisfied and the rest of the TF plane Sq where it is not. 



Given the stationary phase points {(jUj,Tj)}j that are solutions to (17), the stationary phase approximation is 
invoked by replacing ([4]) by: 

x{t) ~ y^ijjx u {T)hl{T-t)e4^r\ (23) 
where S is a subset of the TF plane defined as S = {[J i Si} U So, Si is the neighbourhood of the ith stationary 



phase point such that (/ii,Tj) G Si and [j i Si contains all points in the TF plane that satisfy (19). Equation (23) 



promises sparsity directly from analysis without the computationally expensive re-synthesis step associated with 



most methods for sparse atomic decomposition. The atomic decomposition of (23) is sparse in the sense that 
V(/x, r) 3 S the coefficient X w {t) is implicitly set to zero. However there are no guarantees about the degree of 
sparsity that can be achieved or the quality of reconstruction that might be expected apart from the usual ones that 
might be expected from a well-designed snug or tight frame [20]. This will be explored in the following section. 
Together the analysis steps of Q, ( [TTj ) & ( fT2] ), the selection inequalities ( fT9| ) & ( [20] ) and the synthesis equation 



(23) form what might be called a time-frequency stationary phase approximation (TFSPA). 



V. Elementary Signals 

This Section is devoted to a consideration of the stationary phase regions of the TF plane associated with some 
important elementary signals and also in identifying regions of phase -rate dominance. Because of the elementary 
nature of the signals there is some hope that closed form solutions are possible. These elementary waveforms are: 
an impulse; a single tone; a linear chirp; a decaying phasor. Closed form solutions that define the stationary phase 
points are derived. Within the constraints of the available space results of the numerical evaluation of these regions 
are also provided to confirm and expand upon the theoretical results. 

The results are for order n = 4 gammatone (c = 0) and gammachirp (c = 4) filter banks at the cochlear 
spacing described in Section [if] Numerical evaluation considers a frequency band from Lo m in = 2ir x 100 rad/s to 
Umax = 2tt x 5000 rad/s at a sampling rate of 20 kHz. For the gammatone, = 25 and gfj = § and for the 
gammachirp, g^i = 36 and g£ = ^- These values are commensurate with ERB figures for the cochlea, c.f. HTH . 
The spacing of the filters in frequency is such that there are 4 filters with their maximum gain within the ERB 
of each filter ll20ll . In this case 103 filters are used to cover the stated band. The combined frequency response of 
the analysis and synthesis filters banks has a linear phase and is flat to within a fraction of a dB over the band. 
Thus, over the band of interest, the combination of analysis and synthesis filter banks, ([2]) & ([4]), is effectively 
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distortionless. All niters are approximated by simply truncating them at a suitable point to give finite impulse 
response (FIR) niters. Anticausal niters (e.g. synthesis filters) are simulated by incorporating suitable delays. 

10 and C 2 



For the inequalities ((19 1 and (|20[)), C\ = 10 and C2 = 1. The former was found experimentally to provide a 
good indication of the stationary phase regions for the filter bank described above and for all combinations of filter 
spacing and prototype filter described in Section [ITJ The latter is an extreme value. The inequalities, (20) of Section 
III and (54 1 of Appendix [Bj suggest a larger value such as C2 = 10. However a value of unity was chosen because: 
(i) it makes clear the boundary in the TF plane between regions where the projection p(r, p) is greater than the 
norm 



|a(r, p)\\; (ii) it reflects a desire to test the degree of sparsity that could be obtained. Thus (20) becomes 



Studies of ridgelets and skeletons as in 
possible. 



p(r,p) > ||a(r,/x)|| (24) 
12l suggest that for asymptotic signals, very sparse representations are 



A. An impulse 

The impulse is obviously important because it is an extreme example of a transient signal. Further it is not 
an asymptotic signal in the sense considered in |8] and as such illustrates the advantages of applying to PSP 
to the synthesis integral rather than the analysis integral. For an impulse x{t) = 5(t) the analysis of (|2]) yields 
X u (t) = /3/i(/3r)e^ r . From which, (ITT) and <fl2j( give 



P 



KM 



V 



(25) 



and 




^{2 + /3r^j + Pg(u)r,} + j{/3r - Pg(co)} 



For a Gaussian pulse g{ui) = 0, 77 = and 



KM 

KM 



KM 

Kf3r) 



J_ 

St 



^ j, c.f. Appendix 



A 



(26) 

= —jdr. For a gammachirp pulse g(u) = 5, 77 = — 1/n and 
Thus, for all the three pulse type at filter spacing considered, 



the stationary phase region is defined as r = g(uj)7wui, i.e. a contour in the TF plane at the group delay. 

As to phase-rate dominance, the most straightforward case to consider is a uniform filterbank constructed from 

and a T = [-/3 2 r 0]. Using (20 1 and fi 



t T 



Gaussian filters for which Z T = —(3 2 t, = j^t, 3> t 
the test for phase dominance becomes 

|0 T (r J /i)Wa(r,M)| 



p(t, aO 

Recalling that for uniform filter banks 



du,' 
d/j, 



UJ r , 



|a(r,/x)P 

— Lu m in and thus 



1 duj 



> C 2 



(27) 



P < 



OJr, 



UJ r , 



C-2 



In the following sub-section this result will be considered again and contrasted with a related result involving the 
response to a single phasor. 



More generally it is worth noting that in (26) is a function of normalized time (3r alone and is not dependent 
on the filter frequency 10 since both and ^ are constants for a given filterbank. Z T , on the other hand, is 

dependent on oj indirectly through /3 if the filterbank is non-unform. If Z T varies with uj and Z^ does not then the 
relationship between p{r,p) and ||a(r, p)\\ may also be dependent on u. 

For a gammachirp pulse ( [25] ) and ( [26} become 

1 n + 1 f 1 



Zr 







n 



Pr 



n 



+ jc 



Pr 



(28) 
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and 

respectively. Clearly both 5?{Z T } and ^s{Z^} and hence ||a(r, are not dependent on c. 

Figure fTJa) illustrates the relationship between p(r, u) and ||a(r, fi)\\ for a single filter in a filterbank where 
3w = h H/I \ = ^ anc ^ ^ = rad/s for 3 values of the filter parameter c. The stationary phase point is at 
/3t = n = 4 whereas the peak magnitude of the impulse response occurs is at /3r = n — 1 = 3, c.f. Figure [jjb). 
Insight into the relative behaviour of p(r,fi) and ||a(r, and how that is affected by the value of c, can be 
obtained by considering either a uniform filter bank, i.e. 4^ = 0, or, equivalently, a non-uniform filter bank where 



do; 



<C 1 and <C 1. The latter leads to the approximation 

dw 1 



J i->- 



du/3 



{n - /3r} 



The amplitude rate norm is thus 



l a (^> m)II 



P 



which has a zero at /3t = { 



n-l 
n+1 



n — 1 n+1 

/St n 

}n, c.f. Figure jTj^a). Straightforward manipulation also gives: 

i n r <iu i 



/3c 



/3t 



+ 



d/x/3 



{n - f3r} 



(30) 



(31) 



(32) 



which, not surprisingly, has a zero at the stationary phase point. Consider two limiting cases of (24i. First, as 



/3t —> 0, it reduces to 



> P 



n-l 



{(3t -n} > -{n + l-c} 



(33) 



, which is satisfied if |c| > n — 1. Thus, for example, in Figure jlj a), for 
fir < 3, the inequality is satisfied for c > 4 and phase rate dominance can be achieved to the left of the peak in 
the impulse response. Second, as /St — > oo and c > 0, the inequality reduces to: 

\d»P. 

The LHS is always positive and linearly increasing with /3t, hence a value /3to can always be found above which 
the inequality is satisfied. For {c:0<c<n + l}, increasing c reduces the RHS and hence reduces /3tq. For 
c > n + 1 the inequality is satisfied for all /3t. A degree of caution must be exercised in the application of this result 
as it is based on asymptotic arguments. However it does predict the trend. As c increases from 0, /3to decreases, 
increasing the region {(3t : (3r > /3to} where the inequality is satisfied. This trend is evident in Figure [jja) for 
/3t > 4. In summary, increasing the value of c increases the region to left of the peak response where phase rate 
dominance can be achieved. To the right of the stationary phase point this region can also be increased by increasing 
and positive values of c. The effect of this increase will however, as for the Gaussian pulse, be influenced by the 
relative sizes of ja^} and (3, c.f. (33 1 

Figure [2] shows both the complete TFresponse and results for TFSPA for n = 4 and c = 4. The stationary phased 
points that satisfy (19 1 are indicated in white on Figure |2](b). This stationary phase contour is not co-incident with 
the peak response or ridge but rather lies at the group delay of each filter as shown earlier. Figure |2fb) is obtained 
by only plotting |X w (t)| at points in the TF plane where (24i is not satisfied. As expected from Figure [T| the 
leading edge of the response has been removed. Reconstructions of the input waveform are shown in Figure [3] 
using both full TF plane and TFSPA. In this case they are almost identical despite the removal of the leading edge 
of the response. 



B. A single phasor 

Apply a single tone of A rad/s to the analysis filter bank i.e. x(t) = u(t)e^ xt . In the steady state for t » this 
gives 

X u (r) « e^ XT H(n) 
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Fig. 1. Impulse response of gammachirp filter at 1 kHz: (a) comparison of p(r,fi) and ||a(r, at various values of chirp rate parameter 
c; (b) response. 



where £2 



, with Z T = (3{r/ + and 




{1 + •qPgiu)} 



H(ny L + dJ l * 



(34) 



Thus for all filters considered at the any of the filter spacings considered the stationary phase point occurs when 
A = lo, Vr since, by definition, the group delay of the filter at u is g{uj) — — -g- 
As to phase -rate dominance, the most straightforward case to consider is aga 



s SLauunaiy 

lev f H(0)\ 

2;ain a unifo 



from Gaussian filters. For a Gaussian pulse, g{oS) - 

doj 1 



dfi/3 



o, H(n) 

d(3 



+ n 



and r] 



orm filterbank constructed 
0. Thus Z T = jj3VL and 



Noting that Z T is purely imaginary and Z^ is purely real, ||a(r, //)|| = \Z fJ \ and p(r,/i) = the inequality 

p0|, reduces to 



> 



du 1 



dp 
doj 



(35) 
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(a) 




time (s) 

Fig. 2. Response a gammachirp cochlear-spaced filter bank to single impulse: (a) [X w (t)| in dB; (b) TFSPA - stationary phase points in 
white. 
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Fig. 3. Reconstruction of impulse using full TF plane of Fig[2ja) and TFSPA of Fig^b). 
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For uniformly spaced niters, this will be satisfied if (3 > Note that this in direct contradiction to (27 1 and 

hence it is not possible to design a Gaussian filter bank that exhibits phase-rate dominance in response to both 
impulse and phasor inputs. For non-uniform filter banks the RHS is quadratic in Q. However provided << 1, 

the same general conclusion can be drawn for a large range of and hence A. 

For a gammachirp pulse, g(u) = n/(3, H(tt)/H(Q) = ~ n {^+|^^+ c /"} an( j ^ = _ iy n _ Thus 



1 

n 



and 



c d/3 



It} 



At the stationary phase point where, Q = 0, Z T = —j3/n and Z^ = 0, ||a(r, = (3/n. As |fi| — > oo, 



^ d/j, (5 \ doj \ n do; 



1 



(36) 



As in the previous sub-section, for uniform filter banks or for non-uniform filter banks where ^ <C 1 and ^ < 1 

„ dui n 



d3 



-J 



dp 15 



and hence ||a(r, — > (3/n and 



do; n 
d7^ 



(37) 



(38) 



Figure 4[a) illustrates the relationship between p(r, p) and | |a(r, /u)| | for a single filter in a filterbank where — 



dul 

d M /3 



35 and (3 



l 

du) 9' 

712 rad/s for 3 values of the filter parameter c as in the previous sub-section. As expected 
a(r, p)\\ is approximately constant and virtually independent of c. The 



from the asymptotic analysis, the norm 
projection p{r,fj) is approximately linear in Q, as suggested by (38 1, and much less dependent on c than the 



response to an impulse. For these examples, (24 \ is satisfied and phase rate dominance is achieved outside the 
nominal bandwidth (3 of the filters, i.e. for |Q| > i. The implication is that the output of an analysis filter at ui 
rad/s, in response to a phasor at A rad/s, contributes little to the output of the synthesis filter bank if |A — u\ > % • 
Further, ( [24] ) can be used to locate the vicinity of the peak in the response |X w (r)|. 

Figure [5] shows both the complete TF response and the results for TFSPA. The stationary phase points are 
indicated in white on Figure [5^b). The remaining values of |X w (r)| that are plotted are at points in the TF plane 
where ( 24 1 is not satisfied. The steady state behaviour is well predicted from the theoretical considerations above 
as illustrated in Figure |4] The horizontal white line at 1000 Hz corresponds with the stationary phase point on 
Figure |4] There is only a small region around Q = in Figure [4] where (24 lis not satisfied. This region is just 
visible on Fig. [5] as a horizontal band in red at 1000 Hz. The transient behaviour has similarities to that of the 
impulse response of Figure 2(b). Specifically the stationary phase points, on Figure 2(b) at the the group delay, 
are present at most frequencies on the LHS of Figure [5Jb) apart from a band of frequencies around the applied 
one. Figure [6] shows the reconstruction achieved both with the whole TF plane and with TFSPA. The transient 
response of TFSPA introduces an amplitude modulation that dies away before matching the steady state response. 
This modulation is evident for r < 7 ms, the group delay of the filter, after which the stationary phase points at 
1000 Hz is established. 
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Fig. 4. Steady state response of gammachirp cochlear-spaced filter bank to 1 kHz tone. 



C. Linear chirp 

Consider a linear chirp with a rate 7 rad/s 2 , i.e. x(t) = u(t)e , ^ t ' 1 , where u(t) is the step function at the origin. 
The response is: 



e J 2 



/3/i(/3t)e^V? {r -* }2 (tt 



for r> 4. This integral is intractable and is often approximated using the PSP. However it might appear imprudent 
to proceed with repeated application of that principle. Rather the approach here is to approximate the integral by 
more explicit means. First assume that the frequency jt of the quadratic phase term e'a is approximately constant 
over the temporal extent of h({3t), specifically 7 <C P 2 - Then replace the quadratic phase term with its frequency 



at the filter group delay, i.e. e J 2* « e 779 ^*, to give 
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Fig. 5. TF response of a gammachirp cochlear-spaced filter bank to a 1 kHz tone applied at time zero: (a) |X u (r)| in dB; (b) TFSPA 
stationary phase points in white. 
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Fig. 6. Reconstruction of 1 kHz tone applied at time zero 
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The envelope |X w (r)| will have a peak when — 07 + jt 
points are defined as the solution to Q{Z T } = 0, hence 



H 



H 



-a;+7T-7g(gj) 
— lu+^/t— 79(0;) 




73(07) = 0, thus r = - + 3(07). The stationary phase 

-1 - 73(07) 



P 



and Q{Z^} = 0, hence 



H 



7T + 



-fg(u>) 



H 



-<^+7t— 79(0;) 

3 




o> 



(39) 



Substitution for Q{H(. )/H(.)} gives the stationary phase points as 



07 



+ 



1 + 75(07) 



(40) 



Note that 5(07) is generally negative because 5(07) oc 4 and /3 decreases with 07 (in non-uniform filter banks). 
More precisely, for a Gaussian pulse 3(07) = and hence so is 3(07). For log and cochlear spaced filter banks 
using gammatone and gammachirp pulses g(uj) = 5 and hence 3(07) = — J^gfj, where ^ is positive and ^ < 1. 
Thus, for 7 <C /3 2 , the denominator term 1 + 73(07) lies in the interval [0, 1) and thus the stationary phase points 
lag behind the peak response. When the chirp rate is low, i.e. provided 7 <C /3 2 , the behaviour of the p(r, fi) and 
||a(r, can be inferred from Figures [4] There will be a band in the TF around the stationary phase line where 
the ( |24] > is not satisfied. 

Figure [7] illustrates these points. The applied signal contains a down-chirp (7 = — 10 x 10 4 rad/s/s) which can 
be observed in Figure |7Ja) as a continuous ridge from top left to bottom right. A second lower amplitude up-chirp 
(7 = 8 x 10 4 rad/s/s) produces a ridge from bottom left to upper right. The down chirp produces a series of 
stationary phase points (in white) from top left to bottom right. The position of these points is well predicted by 
(40 1, shown as a black line, apart from bottom right where 7 > 1 . As expected there is a region of the TF plane 



around the stationary phase lines where ( 24 1 in not satisfied. The lower amplitude up-chirp is depicted in a similar 



manner, including stationary phase points, apart from the region where the two chirps cross. In this region the 
larger amplitude chirps hides the trajectory of the lower amplitude chirp. This is a form of simultaneous masking 



that will be explored more thoroughly in Section VI Finally there is a low amplitude artifact that lies between the 
two chirps. This artifact does not contain stationary phase points and could be ignored on that basis in applications 
where signal analysis was the primary objective. The outputs from the synthesis filter banks are shown in Figure [8] 
at the point where the two frequencies cross. TFSPA produces some distortion in the reconstruction when compared 
with that produced using the whole TF plane. 



D. A decaying phasor 

Voiced speech can be viewed as a signal of the form ^ A Ylm Aiii(t — nT)e A '^" nT ' where T is the pitch period 
and {Xi}i are the complex formant frequencies with centre frequencies {^(Aj)}^ and time constants {— Sft(Aj)}j; 
are the complex amplitudes of this atomic decomposition. Because of the relationship to speech it is informative to 
consider the stationary phase points associated with a single atom of this decomposition, i.e. when x(t) = u(t)e xt 
with a single complex frequency A, Jt{} < 0. Given this: 

X u (r) = /3e AT / /i(/3i)e^*e- A *dt 



To proceed further consider the case where h(t) G R and thus restrict consideration to Gaussian and gammatone 
pulses. The normalised time derivative is 



+ A. 
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Fig. 7. TF response of a gammatone cochlear-spaced filter bank to two linear chirps (a) |X^(t)| in dB; (b) TFSPA - stationary phase 
points in white - equation J40l for downchirp in black. 



Note that since h(t) is real, the first term on the right hand side is also real if 3?(A) = u at which point: ^s{Z T } 

h(/3t)e> ut e- xt 



0, 



c.f. ( 1 1 ). In a similar manner 



9 X u (t) 



XJt) 



tdt 
df 



which is purely imaginary for Gaussian and gammatone pulses. Thus setting QjZ^} = for ^{A} = uj gives 

h{/3t)e-^ x)t {t - g(uj)}dt = 0, (41) 



c.f. ( 12 1. While it is unlikely that closed form solutions of ( |4Tj ) for r can be found, nevertheless it is still possible to 
infer something about the nature of the solutions. Because h(f3t)e~^ x ^ t > 0,Vt and {t — g(u)} < 0,Vi < g{oj), no 



dt. 



solution exists for r < g(co). For r > g(u) the integral can be split into two terms, i.e. dt = dt + JL^ ■ 
The first term is a negative constant. The second is is nonnegative monotonically increasing function of r because 
/i(/3i)e _5R ( A )*{i — g(uj)} > 0,Vr > g(u)). Thus, if a solution exists it will be the only solution. Such a solution will 
obviously be dependent on the value of the bandwidth f3 of the analysis filter for which to = 3?(A) and the time 
constant — K(A) of the applied pulse. Thus it may be possible to infer both the frequency and time constant of the 
applied pulse from the position of the stationary phase point. 
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Fig. 8. Reconstruction of chirps in the cross-over region using: (a) full TF plane of Fig. |7Ja); (b) TFSPA of |7jb). 



VI. Simultaneous Masking 



Simultaneous audio masking is the well studied process whereby the presence of one tone prevents the detection 
of a second tone [15]. In this section it is shown that the use of the test (fT9]l leads to a form of simultaneous 



masking. In particular, if the test ( 19 1 is used to indicate the presence or otherwise of a tonal component at uj, the 
amount by which that test is de-sensitized by the presence of a second tone at uj\ is dependent on: (i) the frequency 
separation of the tones {uo\ — uj} ; (ii) the relative amplitude \A\\ of the tones; (iii) the magnitude frequency 

of the filter at uj. Consider a signal x(t) = e ]wt + A\e ]u)lt made up of two tones at uj and 



response 



uj\, where A, the relative amplitude, is positive real. The steady-state response of the analysis filter at uj is 



where H\ 



X u (t) = ^ T + A 1 H 1 e^ T 
H(Qi) and fii = {uj\ — uj}/(3 is the normalized frequency separation with derivative 



lgjj + ii- The time derivative follows: 

d 



dr 



as does the frequency derivative: 



d_ 
duj 



It is convenient to rewrite the time derivative 



and the frequency derivative (14i as 



= -\h ^ + A^H^ 
p duj 

13] ) of the phase of the integrand as 

^{^(t)^(t)}-u;|^(t)| 2 
\Xjr)\ 2 



{X:(r)lXUr)}-g(oj)\XUr)\ : 
\XJtW 



(42) 



(43) 
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/3 2 {\A X H X \ + cos(#)} 2 + {% r jY lj 4Sf { sin W {IAiH^ + cos(e)} + cos(^) sin(#)} 5 



|*(r,/i)|| 2 oc 



{1 + |^i^i| 2 + 2\A 1 H l \cos(9)Y 



(52) 



with denominator 



\XJt) 



and numerator of (42) 



1 + lAji^il 2 + 2|^i^!| cos(^) 







(44) 



where 6 = {oj — uj\}t — 4> and 
one period of the beat frequency {a; — uji}. In (43 1 

%{X*(t)-^X u) (t)} = 



9|X:(r)^-X w (r)|- W |X w (r)| 2 = 

\A 1 H 1 \pn 1 {\A 1 H 1 \ + caa(0)}. (45) 
Z{Ai^i}. The angle variable 6 indicates the position in time with respect to 



(46) 



which, unlike the numerator in (42 1, is dependent on the particular pulse used. 

Consider a gammatone pulse for which g(uS) = —n/f3, H(Q) = and jpS = "^2"^ . Hence 9{i7o} 

-n, 9f{fffffi} = ^gf and fljflb 



-jnHl from which the normalized frequency derivative is obtained using 



(46). 



9{X;(r)|-X w (T)}- 5 H|^(r)| 2 



n 



/31 + fi 2 



\AtHt\ini 

\AtHt 



d/3 



}{|Aiffi|+cos(e)} 



S1i^ + Usin(0) 
da; 



/3i + n? 

After some manipulation this yields equation (52i (see over), where 



sin(^) 



cos(V') 







J-du; 



+ 1 



(47) 



(48) 



(49) 



The numerator of (52 1 is the sum of two oscillatory terms in 6. Assuming that the filterbank is such that (3 3> 
{^f} {hw} I' t * ien ^ e ^ rst osclli atory term, from 9{Z T }, is larger than the second term, from 

To obtain an approximation to the norm, first consider |.Aii2i| <C 1, in which case the denominator of (|52j) is a 



constant. For 3> 1, | sin(*0)| « 1, | cos(-0)| 



du 



(i.e. small and finite) and 



d/3 
do; 



0. Thus the numerator is the 



sum of two oscillatory terms proportional to {|.Ai.ffi| + cos(#)} . Thus the minima of the norm are, approximately, 



at cos(#) 



-A\H\. For 0,i <C 1, |cos(^)| ~ 1 and | sin(-0) | « and hence the two oscillatory terms are 
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approximately out of phase and thus the minima of the numerator can again be found at the zeros of the larger 
term, i.e. cos(0) = —A\Hi. Substitution in (52 1 yields 

II*/ Ml dUJU l^ll \ A l H l\ ,r n , 

ll»(r,rt|U » + m 

For | A 1 S> 1 the first term does not go to zero. Given that this time derivative term much larger than the 
frequency derivative term, an approximation to the norm is obtained by assuming that the latter term is negligible 
and hence 

" lT '^" ~ l + |A 1 F 1 | 2 + 2|^ 1 i? 1 |cos(^)' 
This function of angle has minima at cos(#) = 1, which are: 

ll*(r,M)IUin - +1 (51) 



Together (50) and (51 1 approximate the behaviour of the minima of the norm ||4?(t, as a function of the 
frequency separation Sli and the amplitude of the the second tone as observed at the output to the filter |.Aii2i|. 
When two tones are present the norm will not go to zero. However, despite that, the test (19 1 may be satisfied 
for a given threshold C\ twice per period of the beat frequency at cos(#) = — |Al-HiJ and thus indicate that a 
region of stationary phase is present. It is also evident from ( [50] ) and ( [51) that the behaviour of the norm changes 
significantly at the point where = 1. 

Figure [9] illustrates ||4>(t, /i)|| m in as a function of |^4iiJi| for a particular filter and a frequency separation, Q\ = 1. 
The minimum value of the norm is evaluated both by numerical minimization of ( [52] ) and using the approximations 
of ( [50| ) and ([57]). There is a clear discontinuity at |^4ii?i| = 1 as might be expected from the development of the 
approximations. The approximation are generally a good fit to the minima evaluated numerically. For this example 
the most significant errors can be observed around |.Ail?i| = 1. It is also clear from the numerical evaluation that 
the minimum of the norm is a montotonically increasing function of |-Aii?i| and as such is invertible. Thus if 
a particular value of the threshold C\ in ( [T9] ) is chosen to detect stationary phase points a corresponding value 
of A\H\ can be calculated that will just achieve a minimum of C\. Together ( [50] ) and ( [51) are not guaranteed 
to provide a montonically increasing function. However minor adjustments in the vicinity of the discontinuity can 
overcome this and together they provide a simple means of inverting the function. Figure 10 illustrates how adoption 
of the stationary phase test ( fT9[ ) leads to a simultaneous masking effect. A dB masking tone at 1 kHz is introduced 
that corresponds to the interferer at uj\ rad/s. For each filter frequency u rad/s, (50 1 and ( [5T) are used to calculate 
the amplitude of a tonal component at to rad/s that would produce a value of minimum norm equal to the threshold 
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Fig. 10. Detection threshold of gammatone cochlear-spaced filter bank for C\ = 10 ; dB masking tone at 1 kHz; response of filterbank 
to masker also shown for reference 



value of C\. This amplitude is converted to deciBels to give the detection threshold. For reference, the response 
of each of the filters, i.e. |jt w (r)|, to the masking tone alone is shown. When wi = wa single tone is present and 
the results of sub-section V-B| apply. The norm will go to zero at ui\ and thus the test will be satisfied. There is a 
region around u\ where the detection threshold increases as |u; — uj\\ increases. However outwith that region the 
detection threshold follows the same trends as |JT w (t)| and, as such, this masking effect is not symmetrical and 
affects filters above uji more than below it. 



VII. Conclusions 

The starting point for this paper was an examination of the application of the PSP to non-asymptotic integrals in 
general and TF synthesis in particular. The conclusion was that only one aspect of the PSP, location of stationary 
phase points, is required. The second aspect, approximation of the integral through use of the second derivative of 
the phase of the integrand, is only needed when closed form expressions that approximate the integral are required. 
When this requirement is removed, the second aspect can be replaced with a test for phase rate dominance. Regions 
of the TF plane that pass the test and do not contain stationary phase points contribute little or nothing to the final 
output. Analysis values that lie in these regions can thus be set to zero. In regions of the TF plane that fail the 
test or in the vicinity of stationary phase points, synthesis is performed in the usual way. In re-examining the 
application of the PSP to the TF synthesis integral, a new interpretation of the location parameters associated with 
the synthesis filters leads to: (i) a test for locating stationary phase points in the TF plane; (ii) a test for phase 
rate dominance in that plane. With this formulation the stationary phase regions of several elementary signals have 
been identified theoretically and it has been shown that sparse reconstructions of tones and chirps are possible. An 
analysis of the TF phase rate characteristics for the case of two simultaneous tones predicts and quantifies a form 
of simultaneous masking similar to that which characterizes the auditory system. All of the above was addressed 
from the perspective of deterministic signals. Consideration of random processes is left for future work. 



Appendix A 

This appendix details the relevant characteristics of Gaussian, gammatone and gammachirp prototype filters. All 
have a maximum unit gain at zero frequency. These prototype filters are scaled in time and modulated in frequency to 
form uniformly-spaced, logarithmically-spaced and cochlear-spaced analysis and synthesis filter banks. A Gaussian 
pulse has a non-causal impulse response h(t) = -^==e~* 2 / 2 with a peak at t = 0, a group delay of and a 

derivative = —th(t). It's frequency response is H(£l) = e~^~ with derivative dH d ^ = —QH(£1). An order 
n gammachirp pulse, with chirp rate parameter c has a causal impulse response defined in terms of the complex 
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gamma function T(.) as 



h(t) 



V_^Jji]__ t n- 1 e - { }t +jc ln(i) 



Vt > 0, 



r(n + jc) 

/i(t) = 0, otherwise. This is a minor modification to the gammachirp filter of lfl4l in order to decouple the 
dependency between the location of the peak gain in the frequency response and the chirp rate parameter c. The 
time derivative is: 

dh(t) f n - 1 f 1 1 

= <!-l + ^— +.7C<! t-- 

n 



dl 



h(t) 



The gammatone pulse is a specific case of this when the chirp rate parameter is zero, i.e. c = 0. The gammachirp 
pulse has frequency response: 

hi + i} n+JC 



I n+jc 



The maximum gain of the frequency response is dB which occurs at Jl = 0. The derivative of its frequency 
response is 

dff(fi) {c - jn} 



Mn + *} + i} 



The peak in the magnitude of the impulse response occurs at t 

dZH(0) 

an 



n 



1 while the group delay at Q, = is 



n. The gammachirp pulse, like the Gaussian, is scaled (in time) and modulated to give the analysis 
filter of [3] The chirp rate parameter is unaffected by the time scaling because of the action of the logarithm 
function i.e. jcln(/3i) = jcln(/3) + jcln(t). The time scaling adds a phase shift cm(/3) to the impulse response. 
Thus values of c can be used interchangeably with [14]. From ([52]), the phase derivative of the impulse response 



is: = — f + j, which is zero at the group delay of t = n. Likewise = — ^ at the group delay. 



Consider a complex function /(v) 



vector is V 



f 



9L 9/ IT 



Appendix B 
a(v)e jb ( v } of a vector v 



[v\ V2] T with a,b,v\,V2 G R. The gradient 



Thus 



V/(v) 
/(v) 



a(v) 



" da(v) ' 




' db{v) ' 


dv t 
da(v) 


+ J 


db(v) 


. dv 2 




9v 2 



V a (v)+jV b (v) 



(53) 



The integral of the function around a point vo over an interval S = {v : ||v|| < r} with radius r is Is{^o) = 
Js /( v o + v)dv. Assuming that the function is well approximated by its first order Taylor series over this interval 
and making a change of variable v = Qu, where Q is a orthonormal rotation matrix chosen such that Q V a (vo) — 
[||V o (v )|| 0] T and Q T V 6 (v ) 4 [V w (v ) V fe2 (v )] T , gives: 



/(vo) f {l + u T Q T V a (v )+ju T Q T V b (v )}du 

Js 



/(vo) { J J dnidn 2 



+ {||V a (v )|| +jV W (v )} J j U1&U1&U2 

+ jV 62 (v ) / j niduidt 
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Given that the interval is circular JJ u\du\du2 = JJ U\du2du2, the integral becomes 

s s 



^s(vo) ~ /(v ) < 



duidu2 



+ {l|V„(v )|| +j{V 6 i(v ) + V 62 (v )}} J j uiduid« 2 



By definition, the rotation matrix is given by, Q 



T _ 



cos(0) sin(0) 
- sin(0) cos(6») 



, where 



cos(0) 
sin(0) 



v "( v «) Thus the 
||V.(vo)||- inUS me 



integral over the interval will be dominated by the phase variations if |V&i(vo) + Vf>2(vo)| ~> || V a (vo)||, where 

V 6 i(v ) + V 62 (v ) = [1 l]Q T V b (v ) 

Vl(vo) 



l|V«(v )|| 

and W = [_\ {]. Therefore a test for first order phase rate dominance is 

|Vr(vo)WV & (v 
l|V„(v )|| 



WV fe (v ) 



> ||V a (v ) 



(54) 
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